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Abstract 

The effect of mechanical interactions between cells in the spreading of bac¬ 
terial populations was investigated in one-dimensional space. A continuum- 
mechanics approach, comprising cell migration, proliferation, and exclusion pro¬ 
cesses, was employed to elucidate the dynamics. The consequent nonlinear 
reaction-diffusion-like equation describes the constitution dynamics of a bacte¬ 
rial population. In this model, bacterial cells were treated as rod-like particles 
that interact with each other through hard-core repulsion, which introduces 
the exclusion effect that causes bacterial populations to migrate quickly and 
at high density. The propagation of bacterial density as a traveling wave front 
over extended times was also analysed. The analytical and numerical solutions 
revealed that the front speed was enhanced by the exclusion process, which de¬ 
pended upon the cell-packing fraction. Finally, we qualitatively compared our 
theoretical results with experimental evidence. 

Keywords: Traveling wave, Nonlinear reaction-diffusion model, Bacterial 
colony 


1. Introduction 

In recent decades, much attention has been paid to the collective behavior 
of bacterial populations. This system is used as the prototype for understand¬ 
ing multicellular assemblies, such as tissue and biofilm .lj. The insight into 
the underlying mechanism of dynamics is important to biological and medical 
science. 

To cope with unfavorable environmental conditions, bacterial colonies gener¬ 
ate varieties of pattern formations 0,0. The spatiotemporal pattern formation 
in bacterial colonies results from cell migration and proliferation. These dy¬ 
namics at a continuum level can be described by reaction-diffusion processes 
00 . The simplified model 0 ] relied on a density-dependent (or degenerate) 
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reaction-diffusion equation , which was an extension of the classical Fisher- 


KPP equation 10, 11]. These well-known solutions revealed that bacterial 
density evolves as a sharp traveling wave with constant front speed [2]. In 
our previous work, we found an explicit space-time solution for the generalized 
Fisher-KPP equation in one-dimensional space fl2l ]. This solution evolves from 
a specific initial condition to a self-similar object that converges to the usual 
traveling wave on an extended time scale. Although capable of explaining these 
dynamics, the conventional model omitted the size of the bacterial cell. In real 
systems, most bacterial cells are rod shaped and grow in dense environments. 
Accordingly, the mechanical interactions between cells could play crucial roles 
in the spreading of bacterial colonies. 

Recent experimental and theoretical studies showed that mechanical inter¬ 
actions between cells have important roles in the collective behavior of bacterial 
colonies H- Ts|. The dependence on the elastic modulus of the front speed 
has theoretically been found [19]. It mentions that the migration of bacteria is 


caused by cell pushing rather than self-propulsion in dense colonies [14l . 17, 18|. 


Therefore, we speculate that the exclusion process that prevents the overlap¬ 
ping of cells could play a crucial role in the spreading of bacterial colonies. 
This issue motivates us to extend the conventional density-dependent reaction- 
diffusion equation by incorporating the cell size into the investigation of 
the dynamics of bacterial populations. 

In this work, we considered the systems of bacterial cells growing on a 
thin layer of nutrient-rich fluid medium. The bacteria increased in population 
through cell division and interacted through hard-core repulsion (steric inter¬ 
actions), which resulted in exclusion effects and consequent non-overlapping of 
cells. Although bacteria are self-propelled particles [20||, in colonies of densely 
packed or non-motile cells, bacterial migration was caused by cell pus hing , re- 
sulting from cell growth and division, rather than self-propulsion [14l 1171 Il8|. 
Thus, the bacteria behave similarly as passive particles or nonmotile cells in 


high density environments. Apart from cells, Bruna and Chapman [21] ob¬ 


served that the self-diffusion of hard, spherical Brownian particles in a dilute 
regime decreased as the density increased, due to the diffusion of any single 
particle being impeded by collisions with other particles. However, these colli¬ 
sions encouraged the particle to move toward low-density regions, resulting in 
this biased migration being faster than self-diffusion and enhancing overall col¬ 
lective diffusion. Guided by the work of Bruna and Chapman [Hj , we propose 
that bacterial cells move based on hard-core repulsion and without self-propelled 
motility in dense colonies. 

After incorporating exclusion processes in cell (or particle) dynamics, altered 
diffusion coefficients in the continuum limits were found |22j-|29[. The enhance¬ 
ment or slowing of diffusion depends upon cell length and the available moving 
distance, as shown by lattice-based analysis [28]. In some models, diffusion 
diverges to infinity in closely packed densities [22|, [2^, 0] . Singular diffusion 
has also been modeled through the migration of bacterial biofilm 0EH and 
glioblastoma tumors [32]. However, the effect of diverged diffusion on the prop¬ 
agation speed of cell populations remains unknown. 
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To address this question, we employed a continuum-mechanics approach to 
cell proliferation [33] in order to investigate the spreading of bacterial popula¬ 


tions in the presence of exclusion processes. Additionally, we analytically and 
numerically elucidated the front speed of bacterial colony expansion in terms 
of cell size and discussed the consistency of our theoretical results with the 
experimental evidence. 


2. Continuum mechanical model 

2.1. Constitution equations 

From a macroscopic view, bacterial populations constitute continuum fluid 
capable of reproducing in order to increase cell numbers. By pushing each other 


following cell division l4 [17], 18], population pressure increases as a result 
of collisions between cells and forces cells to move. During movement, cells 
encounter friction from the surrounding fluid medium and the substrate surface. 
For the sake of simplicity, we considered the expansion of bacterial colonies in 
one-dimensional space, regardless of cell orientation. Adapting from [33j| . the 
constitution equations that describe the evolution of the cell density, p(x, t), and 
collective velocity, V(x,t), of the bacterial population at position x and time t 
are given by 


dp 

dt 

— 7 V 


d(pV) 

dx 


+ r(p), 


dp dp dp 
dx dp dx ’ 


(1) 

( 2 ) 


where T(p(x, t)) represents the growth function, p(p(x, t)) represents the internal 
population pressure, and 7 represents the damping constant. Eq. m represents 
the continuity equation with the growth term. We assume that bacterial growth 
obeys the law of population growth as described by a logistic function: T(p) = 
kp( 1 — p/Pm ), where k is the rate constant and p m is the maximum density 
[9j, [33(. Eq. © arises from the force balance between Stokes’ law for friction 
and the pressure gradient, which is similar to Darcy’s law describing fluid flow 
through a porous medium. 

We model the bacterial cells as non-overlapping hard-rod particles of average 
length, a, that interact through hard-core repulsion. In high-density environ¬ 
ments, bacterial self-propulsion can be ignored, since it is dominated by collision 
between cells. This defines the bacterial cell as a passive particle or non-motile 
cell that obeys the laws of thermodynamics. For hard-rod fluid in one dimension, 
the exact pressure is given as 


pk B T 
PKP) = — 


ap 


(3) 


where k B is the Boltzmann constant and T represents the temperature [34l-i36j. 
In our case where bacterial cells behave as passive particles, the temperature 
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relates to the average translational kinetic energy of a cell, < Ek >= (1/2)fcgT, 
we assume that the temperature is constant in our system. The pressure in 
Eq. © diverges to infinity at closely packed density: p —> \/a. Notably, in 
dilute density, p —> 0, Eq. © recovers the pressure of an ideal gas: p = pksT. 
As shown by [37h391 |. the pressure for dilute active particles is similar to the 
ideal gas, except that the source of kinetic energy comes from the swim speed, 
Uq: ksT oc Uq (33,01. As will be shown later, the temperature source is not 
important; as long as it is constant, the dynamics of our model are invariant. 


2.2. Dimensionless equations 

We define the maximum density as p m = l/a m , where a m represents the 
average length occupied by one cell and a m > a > 0. The logistic law limits 
the growth of bacteria, such that 0<p<p m <l/a. For convenience of further 
analysis, we introduce the following dimensionless quantities: 0 < u = p/p m < 
1 , v = [y/( kpmksT )] 1 ^ 2 E, 0 < e = ap m = v/om < 1, t' = at , and x' = 
[(fcy)/ {PmkBT)] 1 ^ 2 ^. In one dimension, the packing fraction, (e), represents the 
length fraction, which is equivalent to the area and volume fractions in two 
and three dimensions, respectively. We then rewrite Eq. © and Eq. © by 
employing Eq. © in dimensionless form: 


du 

~dt 

v 


ox 

1 du 
(1 - euf dx’ 


u), 


(4) 

(5) 


where the prime has been dropped. From Eq. ©, the migration of bacterial 
populations is biased to move down the density gradient and enhanced by the 
exclusion process, implied from the factor 1/(1 — eu) 2 . This factor increases with 
the density and diverges to infinity as e —> 1 at u = 1 , which causes the bacterial 
population to migrate faster at higher density. This singularity has appeared in 
similar models using different approaches [22, 23[ 129H321 140|. Fortunately, the 


velocity in Eq. © is finite, since du/dx —> 0 at u = 1. The density inside of 
the colony reaches a saturated value, except in proximity to the colony edge. In 
this regime, the density distribution is homogeneous and its gradient approaches 
zero. 

Substituting Eq. © into Eq. ©, we obtain a nonlinear partial differential 
equation: 

MO,)®”' 


du 

~dt 


d_ 

dx 


dx 


g(u), 


( 6 ) 


where M(u) = u/{ 1 — eu) 2 and g(u) = u(l — u). Eq. © is in the same form 
as the density-dependent reaction-diffusion equation, however, the migration 
and diffusion coefficients differ. This is unrelated to the mean-square displace¬ 
ment, however, M ~ pdp/dp. In this model, the populations migrate based 
on the collision between cells as opposed to a random walk. A similar coef¬ 
ficient represents the contribution of hard-core repulsion between cells to the 
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migration of myxobacteria in a dense phase [40| ■ Eq. © is degenerate based on 
M(0) = 0, which results in the sharp interface separated between occupied and 
cell-free regions. In a very dilute system (e —> 0), Eq. © recovers the conven¬ 
tional degenerate Fisher-KPP equation [7H9| , for which an explicit solution was 
determined in our previous work [12]. 


3. Traveling-wave solution 

We focused on behaviour of the system over extended times, during which 
the population density propagates as a traveling wave: u(x,t) = 4>(z), where 
z = x — ct and c represent the front speed §]. Substituting the traveling-wave 
solution into Eq. ([6]), we obtain 


t (Wf)+' 

dz \ dz J dz 


(7) 


In the degenerate model, the density must vanish at the finite position, 2 *(< 
oo), that undergoes the sharp interface. We then consider the density pro¬ 
file that satisfies the following conditions: cf>(— oo) = 1 , cj>(z) = 0 for z > z*, 
^</>(—oo) = 0, and A^(j)(z*) ^_0. Additionally, for e £ [0,1), oo)) < oo 

and M(cj)(z )) = 0 for 2 > 2 * [4l|. Multiplying Eq. J7]) by M(</>)d0/ dz and then 

integrating with respect to 2 from —oo to z*, we obtain c M(</>) (jpz'j dz + 




= 0. Under these density profile con¬ 


ditions, the last term on the left-hand side is zero. Finally, we obtain the front 
speed: 


Jo A7(0)g(</ > )# 


( 8 ) 


To obtain the closed-form of the front speed, c, the solution for the density 
gradient, dcf>/dz , is required. 


3.1. Approximate solution 

Although the exact solution of Eq. © remains unknown, we can find the 
approximate solution by employing the perturbation method [42j]. By defining 
w{4>) = d(f>/dz, we rewrite Eq. ([7]): 


M(<j))w ^ + M'{(j))w 2 + cw + = 0, 

d<p 


(9) 


where = dM (</>) / d(j>. The migration coefficient can be written in the ex¬ 

pansion form: M ( (f >) ~ (f> (1 + 2 (f>e + 3 (j) 2 e 2 + •••). We then look for the solution 
of Eq. © in the power series of e: 


ui(</>) = w 0 {c/)) + wi(4>)e + w 2 {(/))e 

= Co + C\€ + C 2 fi 2 + ’ ‘ ‘ , 


c 


5 


( 10 ) 

( 11 ) 





where m;(0) and Cj, that i € {0, 1 , 2 ,..., 00} are coefficients to be determined. 
Substituting Eq. m and Eq. m into Eq. ©, we obtain the equation for each 
order as follows: at e°, 


(fnvo^YT + w o + c 0 w 0 + ^(1 - 0) = 0, 
dcp 


and, at e. 


4>w 0 - 


d(f> 


dcj) 




„ , o dwo 
+20 wq—— 
dcp 


■ 4 0UIq + CiWq = 


( 12 ) 


(13) 


Eq. (TH?1) has the known solutions: wq = (1/a/2)(0 — 1) and co = 1/a/2 0-0,142']. 
Substituting these solutions into Eq. m, we obtain a linear first-order ordinary 
differential equation: 


0(0 — l )~^ + (30- l)^i 

+3v / 20 3 - 5^20 2 + (2^2 + ci)0 - ci = 0. 
After finding the integrating factor |43] , we obtain its solution: 


Wl ( 0 ) = 


1 


C 3 a/2 


( 0 - 1 )' 
ci 7V2 


0 4 + 2v / 20 3 


+ 


(ci + a/2^ 0 


- Cl 


(14) 


(15) 


where C is the integral constant. To prevent the singularity at 0 = 0 and 0 = 1, 
we require that C = 0 and — q. 2^/2 — ^ + (ci + a/ 2) — Ci = 0. 

Thus, we obtain 



(16) 


Substituting Eq. (fl6l) into Eq. m we obtain 


Wl(0) 


- i +(0- 1) W- 1) . 


(17) 


Finally, gathering all terms, we obtain the approximate solutions with the cor¬ 
rection of 0(e 2 ) 


dcj) 



c 


( 1 + r) + ° (e2) - 


(18) 

(19) 
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The density gradient approaches zero when the density reaches the maximum 
value, </>—>■ 1, as expected. By using w(<f) = d(f>/dz, we can calculate the 
approximate density profile: 


</>(-) 


1—exp [b(z—zp)] 
1—a exp [b(z— zq)] ’ 

o, 


z < z 0 
z > Zo, 


( 20 ) 


where a = , b = , and zo represents the initial front position where 

</>(z 0 ) = 0. 


3.2. Front speed 

The front speed is the collective velocity at the edge of the colony, c = 
v(<f>(z*)) = v(0). Based on the correction of 0(e 2 ) from Eq. (fTHl) . the front speed 
increases linearly with packing fraction (e). However, substituting Eq. m into 
Eq. ([8]) and after integration, we can obtain a more precise front speed: 


:(e) = 


5 (4e — 6) In (1 — e) + e 2 — 6e 

V2 e (2e 2 - lie + 8) In (1 - e) - 7e 2 + 8e' 


( 21 ) 


The front speed depends upon the packing fraction of a cell. Therefore, the 
front speed recovers the usual value, that cq = l/y/2 ~ 0.7071, in a very dilute 
regime, as e —> 0 l42| . In a closely packed regime, as e —> 1, the front speed 

approaches a finite value, that c(l) = 10 /a/ 2 ss 7.071, and increases by a factor 
of 10 from the dilute regime. 


4. Numerical results and discussion 

As the correction of our approximate solutions is limited to 0(e 2 ), it is coun¬ 
terintuitive, given that the model is designed for capturing dynamics at high 
density. To obtain the actual results at high density, we solved Eq. © directly 
and subjected the solution to a zero-flux boundary condition using the numeri¬ 
cal method. In Eq. ©, the migration coefficient increases with density, which 
is inefficient when solving with an explicit finite-difference scheme 44]. Unfor¬ 
tunately, solving with the standard implicit-numerical scheme is also difficult 
because of the factor 1/ (1 — eu) 2 . We found that the simplest algorithm that 
overcomes these obstructions is the nonstandard fully implicit finite-difference 
method [3o| . This algorithm has proven stable enough to explore the dynamics 
at high-packing fractions. The detailed algorithm is described in the Appendix. 

Although our model is not expected to be accurate for dilute systems, since 
it has neglected bacterial self-propulsion, we focused on bacterial population 
dynamics by varying the values of the packing fraction, e G [0,1). For our 
computation, we chose the spacing step and the time step, such that Sx = 0.05 
and St = 0.01, respectively. The computations were performed on 3000 grids 
for e € [0,0.5] and on 5000 grids for e G (0.5,0.99], with 8000 iterations. For 
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e = 0.999999, the computation was performed on 120,000 grids with 150,000 
iterations. The initial density profile, uq(x), was set to a step function: 


u 0 (x) 


1 , x <r 0 
0 , x > r 0 , 


( 22 ) 


where ro represents the initial front position. To ensure that it was far enough 
from the boundary at origin, we set vq = 50. 




Figure 1: (Color online) The demonstration of density profiles, u(x, t), evolving 
from t = 0 to t = 80, obtained by using the numerical method. The dashed 
lines represent the initial density profiles. The data are shown for every t = 8. 




Figure 2: (Color online) The front position versus time, corresponding to the 
numerical density profiles in Fig. ©, from t = 0 to t = 80. The data are 
shown for every t = 1. The markers represent numerical values and the solid 
lines represent the fitting lines for the last 50 data points. R 2 is the correlation 
coefficient. 

The demonstration of the density profiles, obtained from the numerical 
method, is shown in Fig. m for dilute systems (e = 0.10) and dense systems 
(e = 0.99). It was observed that the density profile evolved with the sharp 
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traveling wave with unchanged shape. The front position, r/(t), was deter¬ 
mined by the first position where the density fell to zero. Due to numerical 
deviation, we measured the first position where the density was 1 x 10 ~ 6 , or 
u(rf,t) < 1 x 1CT 6 . The front positions were collected for every t = 1. To avoid 
the transient effects of the initial stage, the last 50 data points were selected 
for fitting with the linear equation, rj = ct + tq. The corresponding front po¬ 
sitions of the density profiles in Fig. 0. as a function of time, were fitted well 
using the linear equation, as demonstrated in Fig. ©. This implied that the 
density propagated with constant front speed, which was equal to the slope of 
the linear equation. We checked the accuracy of our algorithm by considering 
the front speed under conditions of e = 0. In this case, the numerical front 
speed was equal to 0.7074, which displayed an error of 0.04% of the exact value 
(c 0 = l/\/2 ~ 0.7071 0153 ). Finally, we explored the dynamics of bacterial 
populations in a closely packed regime. We set e = 0.999999, in order to avoid 
dividing by zero for the factor 1/(1 — eu) 2 when u = 1. In a closely packed 
system, the numerical front speed was equal to 3.8115, which was less than the 
analytically predicted value due to the inaccuracy of the approximate solution. 
The plot of the numerical front speed versus the packing fraction, as compared 
with the analytical curve generated from Eq. (l2lT) . is shown in Fig. ©. We 
found that the front speed increased with the packing fraction and reached a 
finite value as e —> 1. The analytical results agreed with the numerical data 
for the small packing fraction (e <C 1 ), since the correction of our analytical 
solution was only 0(e 2 ). 

Finally, we compared our theoretical results to experimental evidence. From 
the experiments |45l . m, the average (or typical) velocity dependence upon 
the packing fraction of bacterial suspensions was determined. Below a critical 
packing fraction < 1 , the average velocity of bacteria increased with the pack¬ 
ing fraction and reached the maximum value at the critical packing fraction 
43,j43|. Above this critical point, the average velocity decayed to zero as the 


packing fraction approached one, due to the lack of free space. The increased 
front speed relative to the packing fraction observed in our model qualitatively 
agrees with the experimental observations under the former conditions. Their 
observations under the latter conditions were not observed in our results, given 
that the front speed in our model reached the maximum value when the packing 
fraction equalled one, which represents the closest packing fraction for a one- 
dimensional hard-rod system. Nevertheless, our data showed that the numerical 
front speed in a closely packed regime increased by a factor of ~ 5 relative to 
the dilute regime, which qualitatively agrees with experimental observations 
45 . 4f| showing increases in average velocity by a factor of ~ 3 in suspensions 
of spherical-shaped bacteria 0 and in typical velocity by a factor of ~ 5 in 
suspensions of rod-shaped bacteria H 3 ■ 


5. Conclusion 

This study demonstrated the effect of mechanical interactions between cells 
based on the spreading of bacterial populations by employing a continuum- 
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Figure 3: (Color online) The front speed versus the packing fraction, e. The 
dashed line represents the analytical curve generated from Eq. ED and the 
circle markers represent the numerical results. 


mechanics modeling approach. In dense colonies, bacterial migration is dom¬ 
inated by hard-core repulsion between cells, which causes exclusion processes. 
The analytical and numerical results revealed that the expansion speed of bac¬ 
terial colonies was enhanced by the exclusion effect and dependent upon the 
cell-packing fraction. These findings are qualitatively consistent with experi¬ 
mental evidence. 
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Appendix A. Nonstandard fully implicit finite-difference scheme 

We define the discrete density as u" = u(xj,t n ), where Xj = jSx, t n = nSt , 
Sx is a spacing step, St is a time step, j £ {0,1,2 ,..., J}, nf {0,1,2,..., N}, 
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and J and N are integers. We then rewrite Eq. © as 


fr *" +1 

dt 


d ( 9< +1 \ 

— I Mf —^— ] + />” +1 


9a; 


9a; 


(A.l) 


where M" = M(u") = it"/ (l — at") 2 and /" = 1 — it". Using the standard 
discretized scheme for the differential operators, we obtain 


“” +1 ~ + 

St 


1 

<5a; 


Ml 


9 


n+1 


i+1/2 9a;+ +1/2 

A/f ra _^-i/ n + 1 1 4- f n n n + 1 

~ Mj - 1/2 dx j - 1 / 2 ) + U u a • 


(A.2) 


We discretize the remain gradient terms in Eq. (IA.2I) and then have 

1 


»“+' - «y 


M", '-" +1 


n+1 


( fe ) 2 L • J+1/2 v 

M n (n n+1 ii n+1 \ 4- f r Ni n + 1 

M j- 1/2 “ +-i + Jj u j ■ 

The migration coefficient at the mid-grid can be computed by 

1 


M 


j- 1/2 


M +l /2 


(m^! + m; 


(mj* + m. 


3+1) 

2 \ 


(A.3) 

(A.4) 

(A.5) 


Noting that the correction of Eq. (IA.3I) is 0(St , ((5a;) 2 ). After rearranging Eq. (IA.3I) . 
we have 




*X +1 


+ /3>+i =«"• 


(A.6) 


where 


a 


4 — MA/7_ 1/2j 


e? = i-^/; + m(m;_ 1/2 + m; +1/2 ), 


/i = (5x) . 


du 


(A.7) 
= 0 


We impose the zero-flux condition at the boundary grid, saying fl, that n 
Qr = 0. Consequently, u^ = u% +1 and M£_ 1/2 = M^ 1/2 . We 

then rewrite Eq. (1A.6I) . subjected to the zero-flux boundary condition, in matrix 
form; 

A" • U n+1 = U", (A.8) 
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where 


0 


' 9S 2/3q ••• 

„,n r\n on 

a x V 1 p 1 


(A.9) 


0 


a 


n 

J -1 


e 


n 

J -1 


2a} 


1 


0 


n 

J J 


and 

U"=[uff u? ••• w}] T . (A.10) 

According to the boundary condition, 0} = 1 — St/ft + 2/xMA 2 and 9 1 } = 1 — 
Stfj + 2/j,Mj_ 1 , 2 . The numerical density is obtained by solving the matrix 
equation (Eq. (lA.gjl l iteratively. 

To find the stability condition of this numerical scheme, we use von Neumann 
stability analysis: 

u] = (\) n e ikjSx , (A.ll) 


where A represents the amplification factor and k is the wave number [44^ . Sub¬ 
stituting Eq. (1A.11D intoEq. (IA.3I) . we obtain A -1 = 1—Jf/”—/rM” +1 / 2 ( e lkSx — l) + 
/xM"_i /2 (l — e ~ lkSx ), which can be approximated further: 


A » [1 — Stf? + 4 nM™ sin 2 (kSx/2) + O(Sx)] 1 . (A.12) 

A_stable and temporal non-oscillated numerical solution requires that 0 < A < 1 
(3p| . According to 0 < /" < 1 and 0 < M” < oo, without the growth term, 
(/"), this algorithm is unconditionally stable as long as 5x <C 1 [I44j| . With the 
growth term, the solution slowly grows to a finite value as long as St -C 1. Based 
on Eq. (IA. 121) . this algorithm is adequately stable for this type of problem. 
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